Multi-level quantum description of decoherence in superconducting qubits 
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We present a multi-level quantum theory of decoherence for a general circuit realization of a 
superconducting qubit. Using electrical network graph theory, we derive a Hamiltonian for the 
circuit. The dissipative circuit elements (external impedances, shunt resistors) are described using 
the Caldeira-Leggett model. The master equation for the superconducting phases in the Born- 
Markov approximation is derived and brought into the Bloch-Redfield form in order to describe 
multi- level dissipative quantum dynamics of the circuit. The model takes into account leakage 
effects, i.e. transitions from the allowed qubit states to higher excited states of the system. As a 
special case, we truncate the Hilbert space and derive a two-level (Bloch) theory with characteristic 
relaxation (Ti) and decoherence {T2) times. We apply our theory to the class of superconducting 
flux qubits; however, the formalism can be applied for both superconducting flux and charge qubits. 
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I. INTRODUCTION 

Since the famous cat paradox was formulated by 
Schrodinger , the question whether the range of valid- 
ity of quantum mechanics in principle extends to macro- 
scopic objects has been a long-standing open problem. 
While macroscopic quantum tunneling was observed in 
several experiments 13, S 13 j there is less experimen- 
tal evidence for macroscopic quantum coherence. The 
experimental study of macroscopic superconducting cir- 
cuits comprising low-capacitance Josephson junctions as 
a physical implementation of a quantum computer (see 
Ref. 6 for a review) represents a new test for macroscopic 
quantum coherence. On the theory side, the effect of dis- 
sipation on macroscopic quantum tunneling and macro- 
scopic quantum coherence was put into a quantitative 
phenomenological model by Caldeira and Leggett 0- 

The fundamental building block of a quantum com- 
puter 13] is the quantum bit (qubit)-a quantum mechan- 
ical two-state system that can be initialized, controlled, 
coupled to other qubits, and read out at the end of a 
quantum computation. Presently, three prototypes of 
superconducting qubits are studied experimentally. The 
charge {Ec > Ej) and the flux {Ej > Ec) qubits 
are distinguished by their Josephson junctions' relative 
magnitude of charging energy Ec and Josephson energy 
Ej. a third type, the phase qubit 0, operates in the 
same regime as the flux qubit, but it consists of a single 
Josephson junction. In all of these systems, the quantum 
state of the superconducting phase differences across the 
Josephson junctions in the circuit contain the quantum 
information, i.e., the state of the qubit. Since the su- 
perconducting phase is a continuous variable as, e.g., the 
position of a particle, superconducting qubits (two-level 
systems) have to be obtained by truncation of an infinite- 
dimensional Hilbert space. This truncation is only ap- 
proximate for various reasons; (i) because it may not be 
possible to prepare the initial state with perfect fidelity 
in the lowest two states, (ii) because of erroneous tran- 
sitions to higher levels (leakage effects) due to imperfect 
gate operations on the system, and (iii) because of erro- 



neous transitions to higher levels due to the unavoidable 
interaction of the system with the environment. One 
result of the present work is a quantitative estimate of 
the effect of errors of type (iii) by studying the multi- 
level dynamics of a superconducting circuit containing 
dissipative elements. The multilevel dynamics and leak- 
age in superconducting qubits may be related to the ob- 
served limited visibility of coherent oscillations. Previous 
theoretical works on the decoherence of superconducting 
qubits 0, have typically relied on the widely 
used spin-boson model that postulates a purely two-level 
dynamics, therefore neglecting leakage effects. Ref. O 
includes the dynamics of an attached measurement de- 
vice, thus going beyond the standard spin-boson model 
while still making the a priori two-level assumption. 

In this paper, we present a general multi-level quan- 
tum theory of decoherence in macroscopic superconduct- 
ing circuits and apply it to circuits designed to represent 
flux qubits, i.e. in the regime Ej ^ Ec- However, the 
same formalism can be applied to charge qubits. Flux 
qubits have been proposed and studied experimentally 
by several groups El IH El IH El El. The first 
step in our analysis is the derivation of a Lagrangian 
and Hamiltonian from the classical dynamics of a super- 
conducting circuit; the Hamiltonian is then used as the 
basis of our quantum theory of the superconducting cir- 
cuit. While deriving the Lagrangian and Hamiltonian of 
a dissipation-free electrical circuit is-at least in principle- 
rathcr straightforward, different possible representations 
of dissipative elements (such as resistors) can be found 
in the literature. One possibility is the representation 
of resistors as transmission lines [2I |^ [23, i.e. an in- 
finite set of dissipation-free elements (capacitors and in- 
ductors). Here, we use a related but different approach 
following Caldeira and Leggett by modeling each resis- 
tive element by a bath of harmonic oscillators that are 
coupled to the degrees of freedom of the circuit 0,11113 
(see also Refs. l26l I27I for extensive reviews). 

We develop a general method for deriving a Hamil- 
tonian for an electrical circuit containing Josephson 
junctions using network graph theory [28|. A simi- 
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lar approach, combining network graph theory with the 
Caldeira-Leggett model for dissipative elements, was pro- 
posed by Devoret 0| . On a more microscopic level, cir- 
cuit theory was also used in combination with Keldysh 
Green functions in order to obtain the full counting 
statistics of electron transport in mesoscopic systems . 
Here, we give explicit general expressions for the Hamil- 
tonian in terms of the network graph parameters of the 
circuit. We apply our theory to Josephson junction net- 
works that are currently under study as possible candi- 
dates for superconducting realizations of quantum bits. 
By tracing out the degrees of freedom of the dissipative 
elements (e.g., resistors), we derive a generalized master 
equation for the superconducting phases. In the Born- 
Markov approximation, the master equation is cast into 
the particularly useful form of the Bloch-Redfield equa- 
tions .29] . Since we do not start from a spin-boson model, 
we can describe multi-level dynamics and thus leakage, 
i.e. transitions from the allowed qubit states to higher ex- 
cited states of the superconducting system. As a special 
case, we truncate the Hilbert space and derive a two- 
level (Bloch) theory with characteristic relaxation (Ti) 
and decoherence (T2) times. 



II. OVERVIEW AND RESULTS 

Before presenting a formal derivation, we explain the 
main results and show how they can be applied to cal- 
culate the relaxation, decoherence, and leakage times Ti, 
T2, and T]^ of a superconducting qubit. Our theory is 
capable of predicting more than these quantities since it 
can be used to model the evolution of the entire density 
matrix. However, we concentrate on the relaxation, de- 
coherence, and leakage time in order to keep the discus- 
sion simple. For concreteness, we discuss the IBM qubit 
[19] , which is described by the electrical circuit drawn in 
Fig. n The procedure is as follows. 

1. Draw and label a network graph of the supercon- 
ducting circuit, in which each two-terminal ele- 
ment (Josephson junction, capacitor, inductor, ex- 
ternal impedance, current source) is represented 
as a branch connecting two nodes. In Fig. the 
IBM qubit is represented as a network graph, where 
thick lines are used as a shorthand for RC-shunted 
Josephson junctions (see Fig.[2Jl. A convention for 
the direction of all branches has to be chosen-in 
Figs, n a-iid 121 the direction of branches is repre- 
sented by an arrow. 

2. Find a tree of the network graph. A tree of a 
graph is a set of branches connecting all nodes that 
does not contain any loops. Here, we choose the 
tree such that it contains all capacitors, as few 
inductors as possible, and neither resistors (exter- 
nal impedances) nor current sources (see Sec. lIIIBl 
for the conditions under which this choice can be 
made). The tree of Fig. ^ that will be used here is 




FIG. 1: The IBM qubit. This is an example of a network 
graph with 6 nodes and 15 branches. Each thick line repre- 
sents a Josephson element, i.e. three branches in parallel, see 
Figure |5| Thin lines represent simple two-terminal elements, 
such as linear inductors (L, K), external impedances (Z), and 
current sources (Is). 




FIG. 2: A Josephson subgraph (thick line) consists of three 
branches; a Josephson junction (cross), a shunt capacitor (C), 
a shunt resistor (R), and no extra nodes. 



shown in Fig.|3| The branches in the tree are called 
tree branches] all other branches are called chords. 
Each chord is associated with the one unique loop 
that is obtained when adding the chord to the tree. 
The orientation of a loop is determined by the di- 
rection of its defining chord. E.g., the orientation 
of the loop pertaining to Li (large circle in Fig. ^ 
is anti-clockwise in Fig. □ 

Find the loop sub- matrices Fcl, Fez, Fes, 'PkLi 
Fkz, and Fkb- The loop sub-matrices have en- 
tries —1, or 0, and hold the information about 
the important interconnections in the circuit. The 
matrix Fxy determines which tree branches X (ei- 
ther capacitors, X = C, or inductors X = K) are 
present in which loop defined by the chords Y (in- 
ductors, Y — L, external impedances Y ~ Z, or 
current sources, Y — B). In order to find, e.g., the 
loop sub-matrix Fcl for the IBM circuit (Figs. 
and El , we have to identify all loops obtained by 
adding a chord inductor (L). Each column in Fc^ 
corresponds to one such loop. In our example, there 
are two chord inductors Li and L3; the correspond- 
ing loops are the main superconducting loop (large 
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FIG. 3: A tree for the circuit shown in Figure A tree is a 
subgraph containing aU nodes and no loop. Here, we choose 
a tree that contains all capacitors (C), some inductors (K), 
but no current sources (Ib) or external impedances (Z). 



circle) and the control loop (small circle) . Each row 
in FcL stands for one capacitor C; therefore, in our 
example, Fcl is a 3 by 2 matrix. The entries in 
each column of Fcl are 1, —1, or 0, depending on 
whether the corresponding capacitor (row) belongs 
to the corresponding loop (column) with the same 
(—1) or opposite (+1) orientation or does not be- 
long to the loop at all (0). E.g., for our example, 

[cf. Eq. nnni] 



'CL 



The first column says that the capacitor Ci (part 
of Ji) belongs to the large loop (in the opposite di- 
rection, thus -1-1), capacitor C2 (part of J2) belongs 
to the large loop (in the same direction, thus —1), 
while capacitor C3 (part of J3) does not belong to 
the large loop at all. Similarly, the second column 
of Fcl says which of the capacitors are contained 
in the small loop. 

4. Use the inductances (self and mutual) 




Lt 



T 

'-'LK 



^LK 



(1) 



and external impedances Z(u;) to calculate the ma- 
trices Mo, N, m, S using Eqs. (EH) , ijESl , lES) 
and (|66|l : for a single external impedance, also use 
Eqs. H73(l - (|75|l to calculate the function K{u!), the 
coupling strength fi and the unit vector m. The 
block form of the inductance matrix Lt originates 
from the distinction between tree (K) and chord 
(L) inductors; L is the chord inductance matrix (in- 
cluding chord-chord mutual inductances as its off- 
diagonal elements), is the tree inductance ma- 
trix, and Tjlk is the tree-chord mutual inductance 
matrix. The Hamiltonian, Eqs. lf77l) ~ (|5^ . together 
with the bath spectral density J{lu) cx lmK{uj), 
Eq. (|?^ . represents the quantum theory of the sys- 
tem including the dissipative environment. The 



form of this Hamiltonian, in particular the Equa- 
tions H62(l - (|66|l are the first main results of this 
paper. The evolution of the density matrix p of 
the superconducting phases only is determined by 
the Bloch-Redfield equation (|107|l with the Redfield 
tensor given by Eqs. and H114(l . representing 

our second main result. 

5. Find the eigenstates and eigenenergies of the sys- 
tem Hamiltonian Eq. (|78|l and calculate the matrix 
elements of the superconducting phase operators tp. 
In practice, this task is usually done numerically or 
using some approximation. Typically, only a finite 
number of eigenstates is known. 

6. For two given quantum levels |0) and |1), the relax- 
ation time Ti and pure dephasing time can be 
found using Eqs. (|124|l and (|125|) : the decoherence 
time is then given by 
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The leakage rate ^ is given by Eq. H15U|I . 

We have carried out the above program for two cases; 
for the IBM qubitjl95 (Fig. P) in Se c.lVml and for the 

mill (Fi 



Delft qubit [Ij, [Ij (Fig. © in Sec. |TX| For the IBM 
qubit, matrix elements were calculated numerically; the 
relaxation and decoherence times in the case of a current- 
biased circuit are plotted in Fig. 01 For the Delft qubit, 
a semiclassical approach was taken, and earlier results 
by van der Wal et al. ^3 for a symmetric SQUID are 
correctly reproduced. In addition to this, the effect of 
SQUID asymmetries-either in the self inductance or in 
the critical currents of the two junctions-are calculated 
in Sec. IIX CI It turns out that typical sample-to-sample 
fiuctuations of the critical current of about 10% can lead 
to a sizable decoherence rate at zero bias current. 



III. CLASSICAL NETWORK THEORY 

The goal of this section is to derive a classical Hamilto- 
nian for an electrical circuit containing superconducting 
elements, such as Josephson junctions. An electric circuit 
will be represented by an oriented graph Q = (A/", B) , 
see Fig. n for an example. 



A. Graph theory 

An oriented graph 0j Q — (JV, B) consists of N nodes 
TV — {ni, . . . , njv} and B branches B — {61, . . . , bs}- In 
circuit analysis, a branch bt = (?T.a(i) i ?^b(i)) represents a 
two-terminal element (resistor, capacitor, inductor, cur- 
rent or voltage source, etc.), connecting its beginning 
node nQ(i) to its ending node nb(i). The degree of a node 
n £ Af is the number of branches containing n. A loop in 
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5 is a subgraph of ^ in which all nodes have degree 2. The 
number of disjoint connected subgraphs which, taken to- 
gether, make up G, will be denoted P and the subgraphs 
Gi, each having Ni nodes and Bi branches {i = 1, . . . , P), 

where J2i=i Ni = N and = ^- Foi" ^ach con- 

nected subgraph we choose a tree %, i.e. a connected 
subgraph of Gi which contains all its nodes and has no 
loops. Note that % has exactly Ni — 1 branches. The 
Bi — Ni + 1 branches that do not belong to the tree are 
called chords. The tree of the graph G is the union of the 
trees of all its subgraphs, 7^, containing N — P branches. 
A tree of the graph shown in Fig. ^ is shown in Fig. O 
The fundamental loops J^i of a subgraph Gi are defined 
as the set of loops in Gi which contain exactly one chord 
fi G Gi\%- We define the orientation of a fundamental 
loop via the orientation of its chord fi. Each connected 
subgraph Gi has Fi ~ Bi — Ni + 1 fundamental loops, 
i.e. the graph has F = X^ili = B — N + P fundamen- 
tal loops (one for each chord). A cutset of a connected 
graph is a set of a minimum number of branches that, 
when deleted, divides the graph into two separate sub- 
graphs. A fundamental cutset of a graph with respect to 
a tree is a cutset that is made up of one tree branch Ci 
and a unique set of chords. We denote the set of fun- 
damental cutsets of Gi with respect to the tree % with 
Ci. Each connected subgraph has A^^ — 1 fundamental 
cutsets, therefore there are N — P fundamental cutsets 
in total (one for each tree branch). 

We will use two characteristic matrices of the network 
graph, the fundamental loop matrix (i = 1, . . . .F; j = 
1,...,B), 

{1, if hj e Ti (same direction as /i), 
— 1, if 6j e J-i (direction opposite to fi), 
0, iihj(^T^, 

(2) 

and the fundamental cutset matrix (i = 1, . . . , A + P; 
J = 1,...,B), 

{1, if hj G Ci (same direction as q), 
— 1, if 6j e Ci (direction opposite to a), 
0, if6, ^C,. 

. 

By observing that cutsets always intersect loops in as 
many ingoing as outgoing branches, one finds 

fW(f(^))^ = 0, (4) 

By labeling the branches of the graph G such that the 
first N — P branches belong to the tree T, we obtain 

F(^) = (1|F), (5) 

where F is an (A -t- P) x (i? — A^ — P) matrix. Using 
Eq. Q, we find 

F(^) = (-F^ I 1) . (6) 



B. Electric circuits 

The state of an electric circuit described by a net- 
work graph can be defined by the branch currents I = 
(/i , . . . Is ) , where li denotes the electric current flowing 
in branch hi, and the branch voltages V = (T^i, . . . Vb), 
where Vi denotes the voltage drop across the branch hi. 
The sign of li is positive if a positive current flows from 
node ri£t(j) to n^^j) and negative if a positive current flows 
from node j) to ?T.a(i) ! ^ is positive if the electric po- 
tential is higher at node na[i) than at node ^^(i). 

The conservation of electrical current, combined with 
the condition that no charge can be accumulated at a 
node, implies Kirchhoff's current law, 

f(^)I = 0. (7) 

In a lumped circuit, energy conservation implies Kirch- 
hoff's voltage law in the form 

F^-^^V ^ 0. (8) 

External magnetic fluxes $ = ($i, . . . , ^b-n+p) thread- 
ing the loops of the circuit represent a departure from the 
strict lumped circuit model; if they are present, Faraday's 
law requires that 

F(-^)V = *. (9) 

External fluxes have to be distinguished from the fluxes 
associated with lumped circuit elements (e.g., inductors, 
see below). 

We divide the branch currents and voltages into a tree 
and a chord part, 

I = (Itr,Ich), (10) 
V = (Vt„Vch). (11) 

The 2P branch currents and voltages are not indepen- 
dent; the Kirchhoff laws Eqs. iQ and lO together with 
Eqs. (jSj and © yield the following B equations relating 
them, 

Fleh = -Itr, (12) 
F^Vtr = Veh-*. (13) 

As an example, the N + P tree branch voltages Vtr com- 
bined with the B — N ~ P chord currents Ich completely 
describe the state of a network, since all other currents 
and voltages can be obtained from them via Eqs. H12|l and 
(|13|l . However, in the following, we will use a different 
subset of variables, also making use of the B equations 
that are derived from the current-voltage relations of the 
individual branch elements. 



C. Circuits containing superconducting elements 

For the purpose of analyzing electric circuits contain- 
ing Josephson junctions, we adopt the RSJ model for a 
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Josephson junction, i.e. a junction shunted by a capac- 
itor and a resistor, see Fig. |3| We treat the Josephson 
junctions as nonhnear inductors. A (flux controlled) non- 
hnear inductor 28] is a two-terminal circuit element that 
follows a relation between the time-dependent current 
I{t) flowing through it and the voltage V{t) across it of 
the form 



(14) 



where $(i) = V{t) and / is an arbitrary function. For a 
linear inductor, f{x) = x/L, with L the inductance. 

We begin our analysis by choosing a tree containing 
all of the capacitors in the network, no resistors or ex- 
ternal impedances, no current sources, and as few induc- 
tors as possible (in particular, no Josephson junctions). 
We assume here that the network does not contain any 
capacitor-only loops, which is realistic because in practice 
any loop has a nonzero inductance. A network is called 
proper if in addition to this, it is possible to choose a tree 
without any inductors (i.e., if there are no inductor-only 
cutsets) [23 • Again, it can be argued that this is realistic 
since there always are (at least small) capacitances be- 
tween different parts of a network. But we have avoided 
making the latter assumption here because it spares us 
from describing the dynamics of small parasitic capaci- 
tances. We further assume that each Josephson junction 
is shunted by a finite capacitance, so that we are able to 
choose a tree without any Josephson junctions. Finally, 
we assume for simplicity that the circuit does not con- 
tain any voltage sources; however, voltage sources could 
easily be incorporated into our analysis. 

We divide up the tree and chord currents and voltages 
further, according to the various branch types, 

Itr = (Ic,lK), Ich = (Ij,lL,IiJ,Iz,lB), (15) 
Vtr = (Vc,Vk), Veh = {Vj,Vl,Vr,Vz,Vb), (16) 

where the tree current and voltage vectors contain a ca- 
pacitor (C) and tree inductor (K) part, whereas the chord 
current and voltage vectors consist of parts for chord in- 
ductors, both non- linear (J) and linear (L), shunt resis- 
tors (R) and other external impedances (Z), and bias 
current sources (B). Accordingly, we write 



CJ CL * Cfi * CZ -r CB 
KJ KL KR KZ ^ KB 



(17) 



The sub-matrices ¥ xy will be called loop sub-matrices. 
Note that since Josephson junctions are always shunted 
by a capacitor as a tree branch, there are never any tree 
inductors in parallel with a Josephson junction, F^j = 
0. As a consequence, a tree inductor is never in parallel 
with a shunt resistor, Fkr — 0. 

We then formally define the branch charges and fluxes 
{X = C,K,J,L,R,Z,B), 



lx{t) = Qx(t), 
Yx{t) = ^xit). 



(18) 
(19) 



Using the second Josephson relation and Eq. H19|) . we 
identify the formal fluxes associated with the Joseph- 
son junctions as the superconducting phase differences 
if across the junctions. 



$0 



(20) 



where $0 — h/2e is the superconducting flux quantum. 
It will be assumed that at some initial time to (which can 
be taken as to ^ ^00)7 charges and fluxes (including 
the external fluxes) are zero, Qx = 0, ^x = (including 
if = 0), and * = 0. 

The current-voltage relations for the various types of 
branches are 



Ij 


= Icsincp, 


(21) 


Qc 


= CVc, 


(22) 


II 




(23) 


Ik 




(24) 


Vr 


= Rifl, 


(25) 


Vz(c^) 


= Z(lj)Iz(w), 


(26) 



where Eq. H21|l is the flrst Josephson relation for the 
Josephson junctions (flux-controlled non-linear induc- 
tors), where the diagonal matrix Ic contains the crit- 
ical currents Id of the junctions on its diagonal, and 
siny = {sm(pi,sm(p2, ■ ■ ■ jSinifpfj). Eq. (|22|l describes 
the (linear) capacitors (C is the capacitance matrix), 
Eqs. 123|l and l|24|l the linear inductors, see Eqs. H32|l and 
below. The junction shunt resistors are described 
by Eq. (I25II where R is the (diagonal and real) shunt re- 
sistance matrix. The external impedances are described 
by the relation Eq. H2t)|l between the Fourier transforms 
of the current and voltage, where Z(a;) is the impedance 
matrix. The external impedances can also defined in the 
time domain, 

Vz{t)^ f Z{t-T)Iz{T)dT = {Z*Iz){t), (27) 
where the convolution is defined as 



(f *g)(i) = / f (< - T)g(T)dr. 



(28) 



Causality allows the response function to be nonzero only 
for positive times, Z(t) — for t < 0. In frequency 
space, the replacement uj uj -\- ie with e > guarantees 
convergence of the Fourier transform 



Z(t)e*"*di= / Z{t)e"^*dt. 

-00 Jo 



(29) 



In order to obtain Eq. (|23|l for the inductors, we write 

L i^LK 



^LK 



II 
Ik 



= L 



II 
" Ik 



(30) 
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where L and Iik are the self inductances of the chord 
and tree branch inductors, resp., off-diagonal elements 
describing the mutual inductances among chord induc- 
tors and tree inductors separately, and Iilk is the mu- 
tual inductance matrix between tree and chord inductors. 
Since the total inductance matrix is symmetric and pos- 
itive, i.e. v^Ltv > for all real vectors v, its inverse 
exists, and we find 



II 
Ik 



t-1tT f-1 



■'K 



with the definitions 

L = 

= 



L 



LK- 



(31) 



(32) 
(33) 



Note that the matrices L and L^, being diagonal sub- 
matrices of a symmetric and positive matrix, are also 
symmetric and positive and thus their inverses exist. The 
operators L and Lk as defined in Eqs. Ip^ and Ip!^ are 
invertible since L^^ exists. Moreover, since the inverse of 
the total inductance matrix, see Eq. (|31|l . is symmetric 
and positive, its diagonal sub-matrices are symmetric and 
positive, and thus L, L/^- > 0. 



D. Equations of motion 

In order to derive a Lagrangian for an electric circuit, 
we have to single out among the charges and fluxes a 
complete set of unconstrained degrees of freedom, such 
that each assignment of values to those charges and fluxes 
and their first time derivatives represents a possible dy- 
namical state of the system. Using Eqs. (|17H19|) . (|21l - 
I26II , H3U|) , and H31|) , the time evolution of the charges and 
fluxes can be expressed as the following set of first-order 
integro-differential equations 



$0 

27r' 

Qc 



2^^ 



V,7 = fIjC-'Qc, 



(34) 



Ic = -Fcjic sirup - FcflR *fl 
-FcL {h-'^L - L-iLiKL^i*K) 





-FczL^^ * *z - FcsIb, 


(35) 




= VL = F^iC-iQc + F^i*x + *„ 


(36) 




= Vr = Fp^C~-^Qc, 


(37) 




= Vz = F^^C ^Qc + fJ^z^k, 


(38) 




— —^'mFkl^' ^^L + L/fF/fiL ^LlkLj^^ 


*A 




—^'kFkz^'z^ * — LaFxbIs, 


(39) 



where 'Lzi'-^) = 'Z{LL!)/iLU, and where the convolution is 
given by Eq. H28|) . In the equations for the chord vari- 
ables Eqs. ||33I), lEZI), and (jSHl, we have assumed 
that only the loops closed by a chord inductor (L) are 



threaded by an external flux, $ = (0, ^x, 0, 0, 0). In or- 
der to obtain Eq. 139|1 . we have first used Eq. (|30|l . then 
Eqs. itT^ and ((JEJ, and finally Eq. (pTT|l . We can ehminate 
*K by solving Eq. (|^ . 



= -Ik (F/^lL + FkzI^z^ * *z + ^kbIb) 



(40) 



with the definitions 



La- = (1a - La'Fa'lL-ILlaL^') La, (41) 

Fal = Fkl-Lk'-lIk- (42) 

Further knowledge of the structure of F can be de- 
rived from the fact that Josephson junctions are always 
assumed to be RC-shunted, see Fig. |5| If we label 
the tree branches such that the first Nj < Nc capac- 
itances are the ones shunting the Josephson junctions 
(A^C'=iiumber of capacitances, 7Vj=number of Josephson 
junctions) then we find 



Fcj 

Qc 



' CR — 
QCJ 

Qc 



a, 

Nc-Nj 



Nj 



(43) 
(44) 



where C denotes the capacitors which are not paral- 
lel shunts of a Josephson junction. In general, the 
charges of these additional capacitors represent indepen- 
dent degrees of freedom in addition to the shunt capacitor 
charges Qcj = ^oCif/iir. But from this point onward, 
we will study the case where there are no capacitors ex- 
cept the Josephson junction shunt capacitors, Nc = Nj. 
However, the resulting equation of motion H6I|I with the 
definitions Eqs. I|()2|l - H6()() still allows us to describe pure 
capacitors by treating them as Josephson elements with 
zero critical current 1^ and infinite shunt resistance R. 
With this simplification. 



'CJ 



' CR 



= 1, 



(45) 



and the ip and ip can be chosen as the 2Nj generalized 
coordinates and velocities that satisfy the equation of 
motion 



2tt 



(46) 



-IT- ( FclL^^*l + FczL^^ * *z + FcbIb 

where we have used Eqs. (|34() . 135() . and (|40|l . and intro- 
duced Lj^ = 27rIc/$o, and {Y = Z, B) 



FcY 



lL+L-iLiAL^>LAFAL L"\ (47) 



Fey + FclL ^\ilk^jI'LkFky- (48) 



The remaining state variables obey the following linear 
relations, 

LllL"^*l+LlzL^^ = aL(¥'), (49) 



LzlL + LzzL;.^ * *i 



az(¥'), (50) 
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where we have introduced 







(51) 


Lzz 




(52) 


T 


— r KL^Ki^ KZ, 


(53) 


^ZL 


= F]^'zL/fF/fL, 


(54) 


aL(¥') 


= ^F^lS^' + - F^iLifFifsiij, 

ZTT 


(55) 




= TT^F^^cp - F^^^LaFkbIs- 

ZTT 


(56) 



Note that in the absence of dissipation, L^^ 0, 
Eqs. H49|l and H5U|I are holonomic constraints for the vari- 
ables since Eqs. H49|) and H5U|) can be integrated. If 
Lll, Lzz, and 



Li — ^LL ^ ^Lzi^zZ^ZL, 
Lz = ^'ZZ — i^ZL^LL^LZ, 



(57) 
(58) 



are regular matrices, the solution to Eqs. H49|l and H50|) 
is given by 



4>L = LLl'{aL{^)-'LLz'Lzz*^zi<f)), (59) 
*z = Lzlz' * {az{<f) -LzL^llaLiv')) ■ (60) 



Note that in the limit of large external impedances, 
liz 0, the regularity conditions for Lill, "^zz^ L^, 
and \iz all collapse to the condition that Lll be regular. 
The latter always holds in the absence of mutual induc- 
tances between tree and chord inductors, since in this 
case F KL — Fkl and thus L^^ is symmetric and posi- 
tive, so that its inverse exists. Integrating Eqs. H59|) and 
(|60|l from to to t, using the initial condition (all charges 
and fluxes equal to zero), and substituting the solutions 
into Eq. H46|l . we arrive at the classical equation of mo- 
tion for the superconducting phases tp, 



11 27r 27r 

Cip = -Lj^sincp - R"V - Mov? -Md*ifi- -^'^^x - Sis, 



(61) 



r 



with 



Mo 
N 

m 
S 



* CL L7'LL7 F 



-1 

■'L ^^LL^ 

mL^"'"(w)m"^, 

Fez - Fchi^liyi 



'kl^k^kz-, 



= FcB ~Fcl(Lll) F^^L^Fkb- 



(62) 
(63) 
(64) 
(65) 
(66) 



Although the expression H62|l for the matrix Mo is not 
manifestly symmetric, we show in Appendix ^ that it 
is indeed symmetric, i.e. Mj^ — Mq. This property 
of Mo allows us to write the term Moy in the equa- 
tions of motion H61|l as the gradient of a potential, see 
Eq. (|77|l below. The matrices Is/ldiuj) and R contain all 
the dissipative dynamics of ip; if all external impedances 
(shunt resistors) are removed, then L^^ and thus 
Md(w) (R~^ ^0)- A proof of the symmetry of 
the dissipation matrix, M^ — Mj, and a derivation of 
the representation in Eqs. H64() and H65I) can be found in 
Appendix ^ 

Note that the coupling matrix S to an external bias 
current 1b can be obtained from m by replacing Z by 
B. Physically, this means that the external impedances 
Z can be thought of as fluctuating external currents; in 
particular, if a bias current is shunted in parallel to an 
impedance, Fxz = ±Fxb {X — C,K) then we find 
S = ±m. In deriving the equation of motion (|61|) . we 
have assumed that the external magnetic fluxes and bias 
currents become time-independent after they have been 



switched on in the past, ^ 0, Is ^ (t > to). 
In the absence of mutual inductances between the tree 
and chord inductors, La'l ~ 0, Eqs. H62|) - (|66|l become 
somewhat simpler. 



Mo 




(67) 


N 




(68) 


rh 


= Fez -Fcl^i^F'^iIikFkz-, 


(69) 


S 


= FcB ~FciX'1\f'^]^'LkFkb, 


(70) 




= L + fJ^j^LkFkl = Lll- 


(71) 



It should be noted here that from now on, the shunt 
resistors R can be treated as external impedances by set- 
ting M'^ = Md -|- iwR"^; the only reason for treating the 
shunt resistors separately is that more is known about the 
possible arrangement of the shunt resistors in the circuit. 
We will mostly concentrate on external impedances in our 
examples and neglect the shunt resistors, because in our 
examples R ^ Z. If, in turn, the external impedances 
are pure resistors, i.e. Z(a;) is real and frequency inde- 
pendent, then they can be described as corrections to R, 
i.e. R' = R + Z. 

A few important remarks about the form of the ma- 
trix M(j are in order, (i) We know that IsAdit) is real, 
causal (i.e., Md(i) = for t < 0), and symmetric 
Md = MJ (Appendix 0. A dissipative term in the 
equations of motion with these properties can be mod- 
eled using the Caldeira-Leggett formalism 0. (ii) In the 
lowest-order Born approximation, i.e. perturbation the- 
ory in the equation of motion in the small parameters 
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Z~'^ (see below), the contributions to from different 
external impedances are additive, in the sense that one 
can calculate for each impedance Zi separately, while 
Zj^i — *■ cxD, and then add the contributions in order to 
obtain the full coupling Hamiltonian (see Eq. (|82|l be- 
low). In the same manner, the decoherence rates due to 
different impedances will be additive in the lowest-order 
Born approximation. An exact statement (independent 
of the Born approximation) can be made if L^^ can be 
written as a sum in which every term contains only one of 
the impedances Z^, since in this case M^; = 'Yl!i=\ Md.i 
where Nz denotes the number of external impedances 
and 'M.d.iiio) describes the effect of Zi. From now on, 
we will study the case of a single external impedance, 
bearing in mind that in lowest-order perturbation theory 
the results obtained in this way can easily be used to 
describe the dynamics of a system coupled to several ex- 
ternal impedances, (iii) In the case of a single impedance, 
Mrf(w) has the form. 



Md(a') 
m 



= m/v7i = m/|m|. 



(72) 
(73) 
(74) 
(75) 



where i^(t) is a scalar real function, m is the normalized 
vector parallel to m, and ^/]x is the length of the vector 
m (/X is the eigenvalue of the rank 1 matrix mm-^). 

The dissipation free (R, Z ^ oo, Md = 0) part of the 
classical equation of motion Eq. \^\\ can be derived from 
the Lagrangian 



27r 



(76) 



1 27r 
+ -¥p^MoV5+^¥'^(N*, + SIb), (77) 



or, equivalcntly, from the Hamiltonian 

2 



(78) 



where the canonical momenta corresponding to the flux 
variables ^o(p/2tt are the capacitor charges 



IV. CANONICAL QUANTIZATION OF Hs AND 
SYSTEM-BATH MODEL 



In this Section, we quantize the classical theory for a 
superconducting circuit that was derived in the previous 
Section. The conjugate flux and charge variables ip and 



Qc now have to be understood as operators with the 
commutation relations 



$0 , 
Zn 



— ihSj. 



(79) 



In order to include the dissipative dynamics of the classi- 
cal equation of motion, Eq. H61|l . in our quantum descrip- 
tion, we follow Caldeira and Leggett 0, and introduce a 
bath (reservoir) of harmonic oscillators describing the de- 
grees of freedom of the external impedances. We will re- 
strict ourselves to the case of a single external impedance 
coupled to the circuit (this is sufficient to describe the 
general case in the lowest-order Born approximation, see 
Sec. IIII|) . For the Hamiltonian of the circuit including 
the external impedance, we write 



H = TCs+Hb+'Hsb, 
Ti-B 



1 



K 



SB 



E 

Q 



(80) 
(81) 

(82) 



where Hs is the quantized Hamiltonian Eq. (|78|l . de- 
rived in Sec. IIIII Hb is the Hamiltonian describing a 
bath of harmonic oscillators with (fictitious) position and 
momentum operators Xa and Pa with [xcp^g] = ihSaf}, 
masses TOq, and oscillator frequencies Ua- Finally, Hsb 
describes the coupling between the system and bath de- 
grees of freedom, <p and Xa, where Cq is a coupling pa- 
rameter and m is defined in Eqs. (|65|l and (|75|l . The 
term AU{(p) — (m • <^)^ J2a Ca/2r7iQa;^ compensates the 
energy renormalization caused by the system-bath inter- 
action (first term) 0. It ensures that, for a fixed value 

of (f, 



min [Uiif) + HBiixa}) + Hsb(<^, {x„})] = C/(c^), 



or, equivalcntly, that for all ip 



(83) 



min [HBiM) + nsBiv, {xo.})] = 0. (84) 

The term AU{tp) will not be relevant for the Redfield 
theory to be derived below. 

In Eq. H82I) , we have already anticipated the form of the 
system-bath interaction. In order to verify this and to de- 
termine the spectral density of the bath (the masses, fre- 
quencies, and coupling constants will only enter through 
this quantity, see below), we derive the classical equa- 
tions of motion from the Hamiltonian Eq. H8U|I in the 
Fourier representation. The equations of motion for the 
bath variables are 

— Uj'^maXa = —iriaUj'^Xa — CaUl ■ (f. (85) 

Solving for Xa, we obtain 

m • (p 



Xq 



(86) 
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The equation of motion for 9? is 



V. MASTER EQUATION 



dip \$o 
Using Eq. we find 



(87) 



, dU (2^\ cl 



Comparing Eq. H88|) to the Fourier transform of Eq. (|61|l . 
and using the decomposition Eqs. I|72|) we obtain the ex- 
pression 



The spectral density of a bath of harmonic oscillators is 
defined as 



(90) 



combining Eqs. H89(l and H9()|l . we arrive at 



1 /27r 



i^(u;) = - hr - / 



, uj'J{uj') 



(91) 



We now use the replacement K{uo) K{uj + ie), since 
ii'(aj) is a function of the external impedance Z{u}), see 

Eq. isni, 



1 



= lim 



1 



and obtain 



= P- 



UJ — U! 



K{lo) 



1 (271 



(92) 



Comparing the imaginary parts, we have identified the 
spectral function of the bath (up to prefactors) with the 
imaginary part of the function K{lo) derived in Sec. IIIII 
from the theory of electrical circuits, 



J(c^) = -A*(^) ImX(c.). 



(93) 



The real parts of Eq. 192|l agree due to the Kramers- 
Kronig relation for K{lo), 



2 f°° 

ReK{u)) = P / duo 

Jo 



,uj'hnK{oj') 



(94) 



which can be derived from the causality relation K{t < 
0) = 0, following from Eq. if^ . 



Starting from the quantum theory for an electrical cir- 
cuit containing Josephson junctions and dissipative ele- 
ments, Eqs. ((75H52|l . we derive in this Section a general- 
ized master equation for the dynamics of the Josephson 
phases only. The equation of motion for the density ma- 
trix of the whole system (superconducting phases plus 
reservoir modes in the external impedances) is given by 
the Liouville equation, 



(95) 



Following from Eq. (|80|l . the Liouville superoperator £ is 
the sum of the Liouville superoperators corresponding to 
the parts Eqs. {THJl, (jHT|l . and (|5^ of the Hamiltonian, 
C = Cs + Cb + CsB, where Cxp = ['Hx,p\ for X = 
S, B, SB. In order to study the dynamics of the system 
without the bath, we take the partial trace over bath 
modes, 



ps{t)^TrBp{t). 



(96) 



From Eq. (|95|) and with the additional assumption that 
the initial state of the whole system is factorizable into 
a system part ps(0) and an equilibrium bath part. 



Pb = Zg'^ cxp{-I3Hb), 



(97) 



with the bath partition function Zb = Trexp(— /^Ti^), 
/3 = l/ksT being the inverse temperature, we obtain the 
(exact) Nakajima-Zwanzig equation. 



ps{t) = -iCsps{t)~i dt'J:{t^t')ps{t'), {98) 
Jo 

E(t)ps - -iTi-BCsBe-'^^'CsBPs^PB: (99) 



where we have used that the interaction Liouville super- 
operator has the form £sb — Cgg(^£gg where Cgg and 
Cgg are system and bath superoperators, respectively, 
and that TrBi^sBPs) = 0- The projection superopera- 
tors P and Q are defined as 



Pp = (Trsp)( 
Qp = p- Pp. 



' PB , 



(100) 
(101) 



The Nakajima-Zwanzig equation with Eq. 

is a formally exact and closed description of the dy- 
namics of the state of the system ps, but it is rather 
unpractical since it still essentially involves diagonaliz- 
ing the complete problem in order to evaluate the expo- 
nential in Eq. H99|) . However, the problem can be sub- 
stantially simplified in the case of weak coupling, i.e. if 
II-CsbII ^ W^s + ^b\\- We assume that the circuit con- 
tains a finite number of external impedances. As we will 
see below, the weak coupling condition is satisfied here if 



<C 1, and 



UJ 



knT 



« 1, 



(102) 
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hold for transition energies LUij between all possible levels 
i ^ j, where J(w) is given in Eq. H93|) . If the coupling of 
the external impedance is strong, /i w 1, then the con- 
dition l|102l) requires that the involved impedance (resis- 
tance) is large compared to the quantum of resistance, 



(103) 



In the regime of Eq. (|l()2|l . we can expand Eq. (|99|) in 
orders of the system-bath interaction Lsb- Retaining 
only the terms in first order (Born approximation) yields 

^■2{t)ps = -zTrB/:sBe-''3(^s+^«)*/:5i3Ps ® Pb, (104) 

where the projector Q in the exponent can be dropped 
without making any further approximation. 

The master equation Eq. (|98|l in the Born approx- 
imation Eq. 1)104(1 . although much simpler than the 
general Nakajima-Zwanzig equation, is still an integro- 
differential equation that is hard to solve in general. Fur- 
ther simplification is achieved with a Markov approxima- 
tion 

ps{t) - -iCsPs{t)'^^{t)ps{t), (105) 



dt'S2(i') 



(106) 



Markov approximations rely on the assumption that the 
temporal correlations in the bath are short-lived and typ- 
ically lead to exponential decay of the coherence and pop- 
ulation. In some situations, e.g. for 1/fnoise, the Markov 
approximation is not appropriate [6Ll26|. Also, note that 
the Markov approximation is not unique [sof. 

The master equation in the Born-Markov approxima- 
tion, Eqs. H105|l and l|10t)l) . can be cast into the form of 
the Redfield equations 29] by taking matrix elements in 
the eigenbasis |n) of Tis (eigenenergies w„). 



Pnm{t) ^^nmPnmi^^ 



kl 



R-n 



iPklit), (107) 



where p„„ = (n\ps\m), ujnm = - w^, and where we 
have introduced the Redfield tensor, 

Rnmki = / dtTYB{n\[nsB{t), [Hsb(O), \k{t)) {l{t)\p B]]\m) 
Jo 

(108) 

using the interaction Hamiltonian and system eigenstates 
in the interaction picture. 



nsB{t) 



^ g*(Ws+WB)t-^^^g-j(Ws+WB)t 



(109) 



\k{t)) = e''^^\k) ^e'^^'lk). (110) 
Further evaluation of the commutators in Eq. ((108(1 yields 

n.nmkl — ^Im nrrk'^ "nk 2_^^ irrm ^ Imnk Imnk' 

r r 

(111) 

r/™L - dte-'*^-''TTBnsB{t)lm'HsB{0)nkPB, (112) 

Jo 

poo 

rLL = dte-^'^'-^TYB-HsBiOhmnsBitUpB, (113) 



with HsBit) 
using the relation (F*-"'''' 



(n|e**'^«HsBe-**^«|TO). Note that, 

^k^mi ' ^^'2 Redfield tensor 

.(+) 



Imnk / 



can be expressed in terms of, e.g., the complex Fj^^^ ten- 
sor only. For our system-bath interaction Hamiltonian, 
Eq. |jH3), we obtain 



Rcr[+^nfc = (m-¥j),™(m-<^)„fcj(|w„fc|)-: 



ImF 



(+) 



sinh/3|w„fe|/2 ' 
(m • ¥3)im(m • ip)nk X (114) 



x-P duj 







UJ — LOnk COth —— 



VI. TWO-LEVEL APPROXIMATION 

If a system is initially prepared in one of the two low- 
est energy eigenstates (0 and 1) and all rates Rnmki for 
fc, ^ = 0, 1 and n, m 7^ 0, 1 are negligible compared to the 
rates Rnmki for n, m, k,l = 0,1 (a sufficient criterion for 
this being low temperature, f3u>i2 3> 1), then we can re- 
strict our description of the system dynamics to the two 
lowest levels. The 2-by-2 density matrix of the system, 
being Hermitian and having trace equal to 1, can then 
be written in the form of three real variables, the Bloch 
vector 



Poi + Pw 
p = Tr(crp) I i(poi - Pio) 
Poo - Pll 



(115) 



where cr — {<Jx, cry, Cz) is the vector composed of the 
three Pauli matrices. 

By combining the Redfield equation 1(107(1 with 
Eq. 1(115(1 . we obtain the Bloch equation, 



p = u; X p i?p + po, 

with u) — (0, 0, LOoi)'^ , 

^(-^0111 + ^01007 
PO = I -RoiOO + ^0111 

^(-^0000 ^ ^'llll) 



(116) 



(117) 



and the relaxation matrix 



-RoiOl ~^R(] 



R 



•^0101 



0101 
2R'n 



Olio 

R'/ 



pi/ 

-ftfllOl" 



0001 



0110 -^0101 
R'n 



-fi^OllO 
-^0110 



^0100 -^0111 



'-^0100 + -^0111 



■0001 



■0000 ' 



'-R'llll 



where Ri 



'7unkl 



KeRnmki and R'^mki 



IvaRn 



(118) 



rrikl • 



If wqi ^ Rnmki i wc Can make the secular approxima- 
tion, onl y re taining terms Rnmki with n — m = k — I (see 
e.g. Ref.li, 



R'n 



Rf^. 



0101 
nil 

-"•0101 





-n^oioi 

-^0101 





Rn 



0000 



-^1111 



(119) 
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The ofF-diagonal term i?oioi can be absorbed into the sys- 
tem Hamiltonian as a frequency renormaUzation, wqi — 
ujQi — -Rqioi ' ^^cl we are left with the relaxation matrix 



R = 








(120) 



where the relaxation and decoherence times are given by 



1 
1 



1 



-^0000 + -^1111 



— -^0101 — 



2Ti 



= 2Re(r, 
Re(rW 



(+) 

Olio 



-(+) 



1001^ 



(121) 



0000 



p(+) 

^ 1111 



^'^ 00117 



1 

27\ 



1 



— Re(rQQQQ 



^ 1111 



■^^ ooiii- 



Using Eq. (|114|l . we obtain 
1 

Ti 
1 

T, 



^ = 4|(0|m. v|l)|2j(o.oi)coth 
= |(0|m.c^|0)-(l|m.<p|l 



^01 

2kBT' 

2 Ji^) 



2kBT. 

w^O 



(122) 
(123) 

(124) 
(125) 



Typically, can be made to diverge by changing the 
external fluxes until (0|m • y|0) = (l|m • (p|l). It can 
be expected, however, that this divergence will be cut off 
by effects that are beyond the present theory, e.g. other 
noise sources, higher-order corrections, or non-Mar kovian 
effects. 



A. Semiclassical approximation 

Let us assume that the potential C/(y) describes a 
double well with "left" and "right" minima at and 
(p^. Furthermore, for the moment we make a semiclassi- 
cal approximation in which the left and right single-well 
groundstates \L) and \R) centered at tp^ ^ are localized 
orbitals, i.e. they do not overlap each other. Then the 
two lowest eigenstates can approximately be written as 
the symmetric and antisymmetric combinations of \R) 
and 



|0> 



1 

1 

71 



1 



1 



ti'Ol 



,(126) 



t^oi 



t^Ol 



,(127) 



where loqi — ^/ JS? + e^, A = {L\'Hs\R) is the tunneling 
amplitude between the two wells, and e — {L\Hs\L) — 
{R\Hs\R) the asymmetry of the double well. Since \L) 
and \R) are localized orbitals, we can approximate 

{L\^\R)^0, {L\^\L)^^L, {R\^\R)^^^. (128) 
From Eqs. H12fcif) - (|128|l the eigenstate matrix elements are 



(oiv-li) 



1 A 



(0|<^|0)-(l|v|l) « — A^, 

(^01 



(129) 
(130) 



where Acp = (pj^ — cpj^. Finally, the relaxation and pure 
dephasing times for a double-well potential in the semi- 
classical limit becomes 



1 
1 

t1 



— ) |A¥?-m|^ J(woi)coth-^2L., 

t^Ol / ZKb-L 



|Acp-m| 

t^oi/ 



J(c.) 



2kBT. 



(131) 
(132) 



In this semiclassical approximation with localized states, 
the relaxation and decoherence times both diverge if Acp 
can be made orthogonal to m. For a symmetric double 
well (e = 0), ^ oo for ah Aip. 



B. Quantum corrections 

Quantum corrections to the semiclassical approxima- 
tion discussed in Sec. IVI Al can be estimated by taking 
into account the finite spread of the wavefunction about 
its center, using a (approximate) quadratic Hamiltonian 
at the potential minimum 



H=^\QlC~^Qc+i^] <^'^Ln„V), (133) 



where 



L-i = Mo + diag 



COs{ipLM-i) 



(134) 



Rescaling ip and its conjugate momentum Qc 
-vC (fi, 



Qc = \/CQc, 
we obtain the Hamiltonian 



1 2tt 
fi = —=—fi 



a$o 

QCA = \fCiQc\i 



(135) 
(136) 



(137) 



where the inverse LC-matrix is defined as 



n= ( VC "'"l^LiTiVC \ (138) 



and we have diagonalized the 17 matrix, Jl^j = '^l^i- 
The ground-state wavefunction of the harmonic oscillator 
Hamiltonian Eq. I|138|) is a Gaussian centered at the left 
(L) or right (R) potential minimum. 



12 



1/4 



exp 



(139) 



r 



The wavefunction overlap integral between the left and 
right state is found to be 



2h \ 271 



2 Nj 



1=1 



(140) 

Note that in the classical limit, where all capacitances Ci 
are large, the overlap tends to zero, {L\R) 0. Intro- 
ducing the orthogonalized (Wannier) orbitals. 



\L) = 



\R) 



\L)~9\R) 
v/1 - 2Sg + 

\R)-m 

- 2Sg + <?2 ■ 



we can derive the matrix elements, 
{L\ip\R) = 0, 



(141) 
(142) 
(143) 

(144) 



1 — 2Sg + g'^ 



1-2^5 + 5' 



and the difference, 
{L\^\L)-{R\,p\R) 



(147) 

where S is defined in Eq. (|140|l . By replacing \L) and \R) 
by \L) and |^) in Eqs. ((1^ and ((TTTjl . we obtain 

(0|c^|l) « i^(^l + ^^ Ac^, (148) 

(0|c^|0)-(l|v|l) « + (149) 

Note that in this semiclassical approximation using Gaus- 
sian orbitals, both Ti and T^, Eqs. p31|) and p32|) . and 
thus also T2, are renormalized by a factor (1 -I- S"^ /2)^^ ^ 
but for the symmetric double- well (e — 0), is still 
infinite. 



the secular approximation, the total rates for transition 
out of the allowed qubit states |fc) {k = 0, 1) can be writ- 
ten as 



1 



TL,k 



4V|(n|m-<^|A;)|V(u;fc„)coth 



2kBT 



(150) 



As an example, we model leakage by adding two addi- 
tional levels 1 2) and |3) to the allowed logical qubit states 
|0) and |1) and derive the typical rate for transitions from 
|0, 1) to |2, 3) due to the coupling to the environment. In 
analogy to Eqs. H126|l and p27|l . the excited states orig- 
inating from two coupled single- well excited states \L') 
and \R') can be written as 



|3> = 



1 



\R') 



\L')-Jl 



\R') 



^23 



(151) 
(152) 



where CJ23 = VA'2 + e"^ and {L\L') = {R\R') = 0. We 
model the coupling to the lowest two levels by the per- 
turbation Hamiltonian 



H' ^-6i\L){R'\ + \R){L'\ + h.c.). 



(153) 



and denote the energy splitting between the lowest two 
states \L), \R) and the higher energy states \L') and \R') 
with r]. In the regime rj 3> A, (5, e ^ A', the matrix 
elements of the phase coordinate (p in the coupled states 
\h) are found to be 



1^13) 
l<P|2) 



(i|^|2) « i=-./l^Ac^, (154) 
V2 77 V ^ 



V2ri\ ^ 



The dominant leakage occurs with the rate 
2 



1 



|m • Ay)pJ(?7)coth 



V 



2kBT 



(155) 



(156) 



Note that (thermally activated) leakage is not relevant if 
T <C 77, in spite of a finite rate , because the pop- 
ulation of the excited states in thermal equilibrium is 
exponentially suppressed. 



VII. LEAKAGE 



VIII. THE IBM QUBIT 



We can go beyond the two-level approximation, e.g., by 
looking at the leakage out of the two lowest levels. Within 



In this Section, we use the theory developed in 
Sees. IIIIHVII to describe decoherence and relaxation in 
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a superconducting flux qubit design which is currently 
under experimental study by a group at IBM 0. This 
superconducting circuit resembles a dc SQUID, with one 
Josephson junction replaced by another dc SQUID, see 
Fig. ^ The circuit thus comprises three Josephson junc- 
tions in total. This design has the advantage that it 
provides a high level of control. There are three exter- 
nally adjustable parameters; the external magnetic fluxes 
threading the larger (main) loop and the smaller (control) 
loop, and the bias current Ib- 



A. Current biased circuit 

We first study the decoherence due to a current source 
that is attached to the circuit, see Fig. Q It is unavoid- 
able that the external current source will also introduce a 
coupling to an external impedance Z . In our model, this 
impedance is connected in parallel with an ideal current 
source. The impedance Z{ijj) as a function of frequency 
Lo can be determined experimentally 19]. 

We choose the tree shown in Fig. for the graph rep- 
resenting the IBM circuit {N = 6 nodes and B — Ih 
branches) and obtain the following network graph char- 
acteristics (cf. Sec. lIII|l . 



'CL 



' KL 




' cz 



' CB 



-1 1 

-1 



(157) 
(158) 



The linear inductances are given by 





' ) 




L2 




0^ 


is / 




\ 


Li j 



Lkl=0. (159) 



Using Eqs. ^ and ^ with Eqs. (HHZl, ((T^ . ((T^ . 

we obtain the parameters for the Hamiltonian, 



Mn 



N 



-Lg L1+L3 + L4 -L1-L4 |,(160) 



L2 + L2 
— L3 Li + L4 

—L2 — L2 — Li 



(161) 



where 

Bo = L1L2 + L1L3 + L2L3 + L2Li + L3L4,. (162) 

For the dissipative part, we use Eq. and Eqs. lf7H|l - 
(|75|l . with the result 



K{lu) - Bo/B^ 



(163) 



where B^ = L1L2L3 + L1L2L4 + LiL:iL4 + LiL2Lz{lj) + 
LiLsL^iuj) + L2L'iL^{uj) -f L2LiLz{Lo) + Ls,LiL^(uj), 



which allows us to determine the spectral density J(w) 
of the bath using Eq. H93|l . and 

„ _ Ll{Ll + Ll) + {L3Li + L2{L; + Li)f 



2 / L3L4 + L2{L3 + L4) 

m = L1L3 

^oVm I L^L2 



(165) 



Since the bias current in shunted in parallel to the ex- 
ternal impedance, we find S = — m = —y^m. We can 
further simplify the expressions in the case of symmetric 
loops, Li = Li and L3 = L2, 



ALi + L2 



2Ll + L1L2 + 4Lii^ + LaLz ' 



GLj + 4L1L2 + Ll 
{AL1+L2Y 



(166) 
(167) 



m 



-h 4Lii2 + Ll \ 



2Li + L2 

Ll 1 . (168) 



Moreover, if the control loop inductance is much smaller 
than the main loop inductance, Li L2, we obtain the 
asymptotics 



Kiuj) 



1 



(169) 



M « I, (170) 



(171) 




The second approximation for K{lo) is suitable if ^ 
Z{uj)/Li which holds for lu < lolr ~ 150 GHz for Z w 
1000 and Ll w 100 pH. In Fig. 2| the relaxation and 
decoherence times Ti and T2 in this regime are plotted 
as a function of the externally applied magnetic flux <I>c, 
using a numerical solution of Tig. 

There is an intuitive explanation for the simplified re- 
sult Eq. H171|l : both an external bias current Ib and the 
current fiuctuations from the external impedance Z are 
split equally between the right and left half of the main 
loop (the two halves having equal inductances). For this 
splitting of the current (fluctuations), the inductance of 
the control loop is irrelevant, since it is negligible com- 
pared to the inductance of the main loop. The current 
in the left half of the main loop is further split equally 
between the two halves of the control loop (having equal 
inductances). Thus, the ratio of current (fluctuations) 
flowing through each of the Josephson junctions is 2:1:1, 
which is reflected in the coupling vector m for current 
fluctuations from the bath to the superconducting phases 
ip pertaining to the Josephson junctions in the right half 
of the main loop, and the right and left halves of the con- 
trol loop, and also in the vector S describing the coupling 
of an external current to the superconducting phases. 
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FIG. 4: Relaxation time Ti, decoherence time Tb, and pure 
dephasing time for the current-biased IBM qubit as a func- 
tion of the control flux around the point <I>c = 'l'o/2. 
The main flux $ is chosen such that the resulting double- 
well is always symmetric. While diverges at the point 
$c = ^o/'2 where (Olm • v'lO) = (l|m • Ti has a min- 

imum at that point. The inductances for this example are 
Li = L4 = 100 pH (main loop) and 1/2 = 1/3= 4pH (control 
loop). The capacitance and critical current of the junctions 
are C = 0.1 pF and h = 8.5 M (Lj = ^o/2nh = 39 pH). 
The external impedance is assumed to be Z{0) = 2.5 kfi at 
zero frequency and Z{iJoi) = lOkil at the transition fre- 
quency ujoi; the temperature of the external impedance is 
taken to be 30 mK. 



B. Flux biased circuit 

Further control for the system shown in Fig. ^ in addi- 
tion to a current bias line can is achieved by inductively 
changing the magnetic flux through the two loops, see 
Fig. |31 This type of control also potentially introduces 
decoherence due to fluctuations of the external fluxes. 
Another way of looking at this effect would be to say 
that, again, current fluctuations are caused by an ex- 
ternal impedance in the coil producing the flux; sub- 




FIG. 5: The flux-biased IBM qubit. The coil inductance K4, 
can either be coupled to the main loop via a mutual induc- 
tance to Li or to the control loop via Ls. 

sequently, these current fluctuations are transferred to 
the superconducting circuit via a "transformer" , i.e. via 
the mutual inductance between the coil and the super- 
conducting qubit. As in the case of the external bias 
current, the decoherence processes are unavoidable if ex- 
ternal control is to be applied. The method introduced 
above can be used in the same way as before to derive 
the Hamiltonian and the spectral density and form of 
coupling of the dissipative environment. The network 
graph (N ^ 7 nodes, B = 15 branches) shown in Fig. [3 
has the following characteristics, 

FcL = -1 1 , Fez = FcB = 0. (172) 

V -1/ 

Fkl = ^ 0^' ^Kz = -Fkb = ^ _i ^ • 

The structure of the inductance matrix depends on 
whether the external flux is coupled to the main loop 
or the control loop. 

1. Main flux bias 

For an external coil coupled to the main (larger) loop, 
the inductances are 

where Lc denotes the self-inductance of the coil and M 
the mutual inductance between the coil and the main 
loop. 

Since the system without external coupling is the same 
as for the current-biased version, the system Hamilto- 
nian, i.e. the expressions for Mq, N, and Bq are the 
same as for the current-biased circuit, Eqs. H160I) and 
(|151|l . with L4 = Li. The spectral density is obtained 
via Eq. (|^ and the result K{uj) — Bq/B^ where B^ = 
2Li{L2 + L3){L, + L,)- LsM^ + {L3 {L, + L,)) - M^). 
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For fi and m we find 

. 2{Ll + L2LS + Li) 



(L2L3 + 2il(L2+i3))2' 
1 



-(1^2+^3) 

^Ll + Ll + iL2 + L,)^ \ lI 



(175) 
(176) 



We study the following special cases of this result. If 
the control loop is symmetric, L3 = L2, we obtain the 
simpler expressions 

4ii + L2 



K{oj) 



6R'P 



-2 



(4Li+ 2.2)2^ 



m = ^ I 1 

V6 



(177) 



(178) 



1 



If for a symmetric control loop, the control loop induc- 
tance is much smaller than the main loop inductance, i.e. 
for Li ^ L2, we find 

4Li 



K{lo) 



lU! 



ALi{Lc + L,) -2AP Z{uj] 



, (179) 
(180) 



The second approximation for K{u}) is suitable if ti) <C 
Z{lu)/Lc, Z{oj)Li/M'^ . The intuitive explanation for the 
result Eq. 1] 178(1 is essentially the same as above for 
Eq. H171|l . with the difference that the inductively cou- 
pled current fluctuations couple oppositely to the Joseph- 
son junction in the main loop. 



2. Control flux bias 

An external coil coupled to the control (small) loop 
can be described by the inductances 



f2Li 



M/2 



-M/2 l!'^) ' (^^^^ 
(182) 



where Lc is the self-inductance of the coil and M is the 
total mutual inductance between the coil and the control 
loop. Again, the expressions for Mg, N, and Bq are the 
same as for the current-biased circuit. We find K{llj) — 
Bo/B^ withS^ = i-L3M^ + 8Li{L2iLc + L,) + L3{Lc + 
L,) - AP) + L2(4L3(ic + i.) - Af2))/4, and 

2 16i? + Ll- L2L3 + Ll + 4ii(L2 + ^3) 



Ar 



2(^2^3 + 2Li(L2 + i3))2 



1 



(183) 



y/i4Li + L2)2 + (4Li + L3)2 + {L2 - Ls)^ 
—L2 + L3 

-(4L1 + L3) I . (184) 
4ii + L2 



Since the bias current is shunted parallel to the external 
impedance, we find S = — rh = —y/JIm. For a symmetric 
control loop, L3 = L2, we obtain 



2L2 



2L2{L, + L,)-AP Z{u)' 



M = 



AP 
2L|' 



m 



V2 



(185) 



(186) 



The second approximation for K{uj) is suitable if <C 
Zicj)/L„Z{w)L2/AP. 

The result Eq. (|186|l reflects the fact that in the sym- 
metric case, L3 = L2, a control flux bias only affects 
the superconducting phases in the control loop. The two 
phases are affected with the same magnitude of fluctua- 
tions, but with opposite sign. 



IX. THE DELFT QUBIT 

As a further application of our theory, we study deco- 
herence in a superconducting circuit studied experimen- 
tally as a candidate for a superconducting flux qubit in 
Refs.iniTl The circuit consists of a ring similar to a dc 
SQUID but with three junctions, see Fig.|Sl For readout, 
a dc SQUID is inductively coupled to the three-junction 
(qubit) loop. The readout SQUID is current biased in 
order to find its critical current. The value of the critical 
current can then be used to determine the state of the 
qubit loop. 

This circuit network graph characteristics of the Delft 
qubit are (A^ = 8 nodes, B = 20 branches). 



'CL = 



P\ 

-1 





' KL 



* 



1 

1/ 

-1 



' ' 



'cz — 



' KZ 



' CB — 



'KB 





Voy 



, (187) 
1). (188) 




FIG. 6; The graph representation of the Delft qubit. The 
"qubit" loop (right) involves three Josephson junctions Ji 
(i=l,2,3; thick lines) and is inductively coupled to the SQUID 
(or read-out) loop (left) which comprises two junctions and 
can be current biased. 
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We use the following assignment for the inductances 



L Mr 
Mr L 



R 



1 — {^l), i^LK 



(Ml 
\M' 



(189) 

where L and L' = Ll+Lr are the self-inductances of the 
qubit and SQUID loops, respectively, and M = Ml+Mr 
is the mutual inductance between the two loops. The 
self-inductance of the SQUID loop and the mutual in- 
ductance are divided into parts Ll and AIl correspond- 
ing to the left half of the SQUID loop and parts Lr 
and Mr corresponding to the right half of the SQUID 
loop. We introduce the following notations and con- 
ventions. The Josephson inductances of the five junc- 
tions are given by Lji = Lj^2 = Lj and Lj^j = Lj/(3 
for the three qubit junctions, and Lj-l^r — L'j for 
the two SQUID junctions. The superconducting phase 
differences across the five junctions are denoted with 
ip = {ipi, ip2, V'i, 'PL, Vr), and the capacitances of the five 
junctions are C = diag(C, C, C, C", C"). The externally 
applied fluxes threading the qubit and SQUID loops are 
described by the vector = {^x, ^'x)- In the symmetric 
case, Ll = Lr = L'/2, Ml = Mr = M/2, we obtain 



Mo 



1 



LL' - Af 2 



N = 



1 



LL' - Af 2 



L' 
L' 
L' 

-Af 
M 

-L' 
-L' 
-L' 
M 

-M 



L' 
L' 
L' 

-M 
M 

M \ 
M ' 

M 
-L 
Lj 



L' 
L' 
L' 

-M 
M 



-M 
-M 
-M 
L 
-L 



M^ 
M * 

M 
-L 
Lj 



,(190) 



(191) 



for the Hamiltonian and 



Z{uo)+iuoL'/A' 



(192) 



M = 1/2, m = (0,0,0, l,l)/\/2, and S = -^m. 

Instead of quantizing the classical Hamiltonian 
Eq. JHOll, with Eqs. (|190(l and (|191|l we wiU hnearize the 
dc SQUID in order to separate the degrees of freedom 
that become very massive under the influence of the ex- 
ternal impedance from the other, light degrees of free- 
dom. Subsequently, we will only quantize the light de- 
grees of freedom, viewing the massive degrees of freedom 
as part of the environment. 



A. Linearization of the dc-SQUID 

We start by linearizing the uncoupled (Af = 0) SQUID. 
The equations of motion for the SQUID are 



C 



{•PL - Vr) 



1 

2L 



- (sin (pL - sin pr) 



1 

T' 



- VR- 27r 



(193) 



{VL + Vr) 



2L 



- [sm(pL + simpR) + —If. 



Now we make the expansion 



(194) 



(195) 



where Pl.r denotes the steady-state solution of the clas- 
sical equations of motion Eqs. (|193|) and (|194|l . We first 
find this steady-state solution in the absence of a bias 
current, Ir = 0, using L' <C L'j and assume $fp ^ $o/2, 
with the result 



-(0) 



-(0) 

-'Pr 



- X 



(196) 



Next, we allow a finite but small bias current Ir ^ I'c 
^o/2TrL'j, and with ipL.R. ~ + ^^l,r. we find 



L'jIb 



$ocos(7r$^/$o)' 



(197) 



Starting from the steady-state solution, we can now de- 
rive the linearized SQUID dynamics d(p{t). We assume 
that the external impedance Z{lu) contains a sizable 
shunt capacitance Csh ^ C" and that w <C 1/ \JL'C' 
(« 1500 GHz for typical values C" = IfF, L' = 10 pH). 
Under these assumptions, the effect of the external 
impedance Z « l/iwCsh is to make the coordinate 
+ PR very "massive", i.e.. 



iijjCsh 



K{t~t')p+{t') 



(198) 

the "mass" being C + Csh/4 ~ Csh/4. In order to elim- 
inate PL + PR from the classical equations of motion, 
we introduce p± = ± pR and expand Eqs. H193() and 
(|194|l about the steady-state solution, p± = ip>± + Sp±, 



2L', 




2T7- ^+iiK{^)^ 



6p+{^)) 

(199) 

where we have used the steady-state solution to define 



C+ = cos PL -\- cos PR = 2 cos TT 



COSpL — COSpR, 



-27r 



$0 

L'jIb 



(200) 

$' 

tan TT^. (201) 

$0 



Neglecting C" <C Csh in the equation of motion for (5<^+, 
we can solve for 5p+ (neglecting higher powers of Ir), 



5p+{uj) = -c- c+ + 



iijjL'j/2 



Z{uj) + iLuL'/A 



Sp^ (uj). 

(202) 
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Substituting this back into Eq. (|199|l . we obtain the fol- 
lowing damping term in the equations of motion for b^p^ 



2L 



Z{u:) 



(203) 



with the effective SQUID inductance Lj = 
Lj/4cos(7r$j;/$o) = L'j/2c+ and the effective ex- 
ternal impedance 



— tanTT-^ 



(204) 



where is the critical current of the SQUID junctions 
and the total impedance (heavy SQUID degree of free- 
dom in parallel with external impedance Z) is defined 
through 



1 



iiuLj Z{uj) + iujL'/A 



(205) 



The effective external impedance Z is much larger than 
ujL'j for Ib < I'c or for sin7r<I>y<i>o w 0. Thus, unlike 
the remaining degrees of freedom (including 5lp_) 
are weakly affected by the effective external impedance 
and will be described as quantum mechanical degrees of 
freedom. 



B. Description of the light degrees of freedom 

After having eliminated one degree of freedom from 
the SQUID, the remaining four degrees of freedom 
ip = {Lpi,ip2,^i,5ip-) will are now described by the 
Hamiltonian Eq. I|80|l with the capacitances C = 
diag(C, C, C, C"/2), the Josephson effective inductances 
= diag(L7\L7\/3L7\0), and 



Mo 



1 



LU - AP 



v 


V 


L' 


-M 


L' 


L' 


L' 


-M 


L' 


L' 


L' 


-M 


M 


-M 


-M 


L 



1 

Tj 



0000 
0000 
0000 
0001 



N 




(206) 



(207) 



Since the part of the circuit that was coupled to the bias 
current is described by Z , there is no coupling to a bias 
current left, S = 0. By inspecting Eq. H203|l . we find 
m — (0, 0, 0, 1), ;U = 1, and 




FIG. 7: Circuit representation of the Delft qubit with a lin- 
earized SQUID. The analysis is simplified due to the absence 
of tree inductors (K). 



The results Eqs. H206|) - (|208|l for the reduced system can 
also be obtained from the circuit drawn in Fig. □ with 
C4 = C"/2 and the inductance matrix 



L M 
M L' 
Lj 



(209) 



Using Eqs. fTM|l . (fH^ . (UlI), and we obtain 



1 



1 



LUQl I \2tT J UJQI 



X tauTT^^ coth , 



e 

/27r/B 



£0 
27r 



|m- A<^|2ReZt(0) x 



, tauTT— 



2kBT. 



(210) 



(211) 



We make the approximation that the SQUID is com- 
pletely classical; solving the classical equation of motion 
Eq. Ht)l|) for 6lp- using Eq. H206|) . we obtain the station- 
ary classical solution for Sip^ (with L', L <C L'j, Lj), 



^V-=-^{MI + K)., 



(212) 



where we have used that J^ifi — — 27r<I>/$0: where $ 
is the flux threading the qubit loop, and ^ — = LI, 
where / is the current circulating in the qubit loop. The 
difference between the two minima <p (localized states |0) 
and |1)) is then is constant) 



27r 

m • A(/j = A6lp- w -—M{Il - Ir), 



and since 1^ 



-Ir = /, 
m • Aip : 



An. 

$0 



(213) 



(214) 



Substituting the above result into Eqs. H21()|l and H211|l . 
we obtain 
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X HeZt (wqi ) coth 



2kBT 



(215) 



1 



,^01 

xReZt(0)2A:BT, 



%2 1 f^TTMI 

/„ tanTT— I X 



(216) 



which agrees with earHer results We also obtain an 
estimate for the leakage rate from Eq. p56|) . 



1 

















c. 



lf2TiMI 

tan TT— 1 ReZt(7]) coth 

'o 'I'o 



2kBT 
(217) 



C. Asymmetric SQUID 

Up to now we have assumed that the SQUID ring in the 
Delft qubit is symmetric in two senses; namely that both 
the self- inductances of the left and right halves of the ring 
are identical and the critical currents of the Josephson 
junctions in the left and right halves of the SQUID ring 
are identical. Both symmetries are certainly broken to 
some degree in real physical systems. Below, we study 
both cases, i.e. the case where the self-inductances of the 
left and right halves of the ring are different (geometrical 
asymmetry) and the case where the two critical currents 
are different (Jimction asymmetry). 



1. Geometric asymmetry 

We analyze the Delft qubit again with the inductance 
matrix Eq. H189|l and the asymmetric inductances 



Ml = 
Ll = 



a\ M 
2)t' 



1 



L, 



1 



Y' 



where Ml + Mr ^ M and Ll + Ln = L'. 
the SQUID, we obtain the result 



(219) 
By linearizing 



— tauTT — 

I'c $0 



a — smTT — 

I'c *0 



(220) 

This result implies that if a 3> Ib/I'c^ the decoherence 
rates scale as als/Ic instead of {Ib/I'c)"^ for a <^ Ib/I'c- 
Therefore, for very asymmetric loops, a ^ Ib/I'c mod- 
erate bias currents can already cause large decoherence 
effects. 



2. Junction asymmetry 

For asymmetric critical currents, or, equivalently, 
asymmetric effective Josephson inductances. 



ij;L,fl = ij(l±7/2), 



(221) 



we repeat the linearization of the SQUID keeping contri- 
butions of lowest order in 7. Setting Ib ~ 0, we obtain 



1 



'J-.LM. 



1 
17 



It 



0(7^). 



(222) 



The steady state of the SQUID is then determined by 
the following equations, 



'PL 

1 - ?1 sini^L 
We make the ansatz 



(l-i)sin^. + (l + |) 



1 



$0 

sin (fR 



^ X 



(223) 
(224) 

(225) 



and obtain the result q = — tan(7r<I>^/<I>o)/2; and finally, 



7 K 
- tauTT— . 

2 $0 



(226) 



Comparing this to Eq. p97|l . we see that in order to 
obtain the Redfield tensor and the decoherence rates at 
zero bias current in the presence of a junction asymmetry 
7, we simply have to make the substitution 



7 sm vr 



X 

$0' 



(227) 



Typical values for the junction asymmetry due to pro- 
cessing inaccuracies are fairly large, 7 ~ 10%. The ef- 
fect of a junction asymmetry is more severe than the ef- 
fect of a geometrical asymmetry because for asymmetric 
junctions, decoherence occurs even for zero bias current 
/s = 0. 
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APPENDIX A: SYMMETRY OF M 







In this Appendix, we prove that the matrix Mq de- 
fined in Eq. H62() is always symmetric, Mj^ = Mq. This 
property is required in order to find a potential U{ip) gen- 
erating the force term — Mo<p in the equation of motion. 



We write Mq = Ycl^^^^^cl with 



V = L^^L = 1l 



XLr-F 



W 



(Al) 
(A2) 
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with the ofF-diagonal block of Lt from Eq. (|31|l . 

X = -L-'LLKh],\ (A3) 
= -L^iLl^L-i, (A4) 

and show that VW^^ is symmetric, thus proving that 
Mq is symmetric. Note that in Eq. (|A4|) we have used 
that Lt is symmetric. As a first step of our proof, we 
note that the symmetry of VW~^, 

(VW-i)'^ EE (W^)-iV'^ = VW-\ (A5) 

is equivalent to the relation 

V^W = W^V. (A6) 

As a second step, we use Eq. H42|l to show 

W - W^V = F^iL^ (Y - Y^) LkFkl , (A7) 

I 



where Y = (L]^)~^ — F/f^X. The third and last step of 
the proof is to show that Y is symmetric, i.e. Y = Y-'". 
For this, we rewrite Eqs. (|41|l and (|42|l as 

Ik = (lK-L;^rKLX)-iLif, (A8) 
Fkl = F;^L-X^L. (A9) 

Using these relations, we can show that Y — Lj} + 
X^LX which is manifestly symmetric. This concludes 
the proof that Mq is symmetric. 



APPENDIX B: SYMMETRY OF Md 

From Eqs. and (EOJ), we obtain hEq. 

with 



Ma{u;) = FcL (Ii'LLiILlz - L-ILlkL^'L^^F^z) tz\^)'i^zLl^llF\ 



T 
CL 



(Bl) 



where we have used the identity hj^^ {Lj)'LLzl'zzi^) ^ 
L^]^Li2L^^(ti;). This expression is a quadratic form in 
Fcl and Fez, 

(B2) 

with the definitions 

A = — L^^LL^^Ll^ + L^^LlaL^"^La:Fa-2, (B3) 
B = -Lzl'LII = -F^-^La'FalL^].. (B4) 

Next, we show that must be symmetric, Mj — M^, 
and therefore = Lz and A = B-^. 

The argument for the symmetry of is as follows. 
We consider a generalized model in which the exter- 
nal impedances Z and the linear inductances L and 
K are treated on an equal footing. For this purpose, 
we allow mutual impedances (generalized mutual induc- 
tances) between Z and K and include Z into L by al- 
lowing frequency dependent linear inductances and writ- 
ing Lzii^) ~ Z{uj)/iuj. This leaves us with the follow- 



ing types of circuit elements; tree elements are either 
capacitors C or linear impedances K where Lif(<^) = 
'LK{'^)/i^, branch elements are Josephson junctions 
(non-linear inductors) J, linear impedances L where 
L(cj) = 'LL{Lo)/ioj, and external bias currents B. In addi- 
tion to this, there can be frequency-dependent linear mu- 
tual impedances 7ilk{^), where Lla'(w) = 7jLK{'^)/i^, 
between the L and K branches. The equation of mo- 
tion Ht)l|) can now be derived exactly as before, but in 
the frequency domain, the result being Eq. (|61|l without 
the 'M.d{u}) term, since there are no Z branches. These 
new equations include dissipation which is described by 
the (now frequency-dependent) MQ(a>), the prime distin- 
guishing it from the "ordinary" Mq (see above). The 
matrix Mf, is formally identical to Mg, up to frequency 
dependencies which are irrelevant for the symmetry of the 
matrix. We have shown in Appendix 1X1 that = Mq; 
this proof also goes through for Mq, thus Mq"^ — Mq. 
Since M.'q{oj) = Mq + Md(cj) and both Mq and Mf, are 
symmetric, we conclude that also Mj = M^. Introduc- 
ing m = Fcl + FczA = Fcl + FczB^, we can now 
write Md in the form given in Eqs. Ht)4|) and H65|l . 
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